The stochastic matching problem 
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The matching problem plays a basic role in combinatorial optimization and in statistical mechan- 
ics. In its stochastic variants, optimization decisions have to be taken given only some probabilistic 
information about the instance. While the deterministic case can be solved in polynomial time, 
stochastic variants are worst-case intractable. We propose an efficient method to solve stochastic 
matching problems which combines some features of the survey propagation equations and of the 
cavity method. We test it on random bipartite graphs, for which we analyze the phase diagram and 
compare the results with exact bounds. Our approach is shown numerically to be effective on the 
full range of parameters, and to outperform state-of-the-art methods. Finally we discuss how the 
method can be generalized to other problems of optimization under uncertainty. 



One important aspect of the statistical physics ap- 
proaches to disordered systems is the broad range of their 
interdisciplinary applications. Systems with frustration, 
structural disorder and uncertainties are in fact ubiqui- 
tous in many fields of science and their study has greatly 
benefited from the algorithms which have emerged at the 
interface between statistical physics of disordered sys- 
tems and computer science. 

One of the key problems has been the so called match- 
ing problem [1], which for the case of random instances 
was among the first to be solved by statistical physics 
methods (2j and later by rigorous mathematical tech- 
niques [3]. Matching is a constituent part of many prob- 
lems in different fields, ranging from physics (dimer mod- 
els [4]), to computer science (vision [5]), economics (auc- 
tions [B]) and computational biology (pattern matching 
[7]). It can be formulated simply (given a graph, find 
the largest possible subset of edges without common ver- 
tices), and is of polynomial complexity [S]. 

The stochastic version of matching is a basic example 
of optimization under uncertainty [HI llOj , which consists 
in finding the minimum of a cost function depending on 
some stochastic parameters, given just some partial infor- 
mation about their value. Most real-world optimization 
problems involve uncertainty: the precise value of some of 
the parameters is often unknown, either because they are 
measured with insufficient accuracy, or because they are 
stochastic in nature and determined only after some deci- 
sions have been taken. The objective of the optimization 
process is thus to find solutions which are optimal in some 
probabilistic sense, a fact which introduces fundamen- 
tal conceptual and computational challenges. Stochastic 
matching problems are in fact known to belong to higher 
computational complexity classes ranging from NP-hard 
to PSPACE-complete [S] depending on how stochasticity 
is introduced. 

Here we apply a new method for stochastic optimiza- 
tion problems to the two-stage matching problem. This 
new method, which builds on the formalism of Survey 



Propagation (SP) jTTHT3] and of the cavity method, is 
partly analytic and allows to optimize the expectation 
of a stochastic cost function by estimating the statistics 
of its minima, without resorting to explicit (and costly) 
sampling techniques. 

In the following we define the problem, describe the 
method we propose for solving it, and discuss its phase 
diagram. We find that for large connectivity the problem 
enters a computationally "hard" phase where standard 
heuristics fail. In particular, we perform a detailed com- 
parison with Stochastic Programming using state-of-the- 
art solvers. While our method has a good performance in 
both phases. Stochastic Programming turns out to be im- 
practically slow in the region of large connectivity, and to 
have a significantly worse performance than our method 
in the region of small connectivity. Finally, we report 
about applications to problems which are NP-hard also 
in their deterministic setting. 

The two-stage stochastic matching problem. We 
study a variant of the problem introduced in |14H16j . 
where it is shown to be NP-complete. We are given a 
bipartite graph G — {L, R, E) with L further partitioned 
in Li and L2, and a set of independent probabilities 
P = {Ph G ]0il[j'2 S L2}- The nodes in Li are de- 
terministic, while the nodes in L2 are stochastic: I2 G L2 
will be available for matching with probability pi^ . In the 
first stage the nodes in Li are matched, knowing only the 
probabilities p. In the second stage, the available nodes 
in L2 are extracted according to p and they are matched. 
The objective is to maximise the size of the final match- 
ing. 

We introduce two sets of binary variables, Xi = 
{xi^r e {0,1}, (^ir) e E : li e Li} and X2 = {xi^r & 
{0,1}, (^2'') & E : I2 ^ L2}, to represent the possible 
M C E, with xir = 1 iff (Ir) e M. We also introduce a 
set of binary parameters t — {t;^ S {0, 1}, ?2 G L2} with 
ti^ = 1 iff ^2 is available for matching in the second stage. 
We define an energy function £(xi,X2,t) counting the 
number of unmatched vertices among the available ones. 
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The problem consists in finding 

= argminlEt min£(xi, X2, t) (1) 

XI X2 

subject to the matching constraints J2iedr ^ir < 1 (Vr £ 

■S)' J2redh^hr < 1 (V/l e Li) and Erea/^ 2^'2r < 

t2 (VZ2 G -^2), where dr = {I e L : (h) e -E} and simi- 
larly for dli and dl2- 

The main difficulty of the problem stems from the fact 
that Et minx2 f (xi, X2, t) has a highly non trivial depen- 
dence on Xi. In order to overcome this difficulty, we 
shall use the cavity method to first compute the mini- 
mum energy relative to X2 for fixed xi and t, and then 
to compute the average over t of this quantity. 

Minimizing relative to X2 for fixed Xi and t. Once 
Xi is determined and the stochastic parameters t are 
set, it is straightforward to find the optimal X2. A 
possible way of doing this is by Max-Sum (MS), as 
discussed in [T7]. We introduce the cavity fields ui^r 
propagating from I2 & L2 to r G R and hi^r prop- 
agating in the opposite direction. The MS equations 
are ui^r = — max[— 1, maxr'gai2\r ^bi"] and hi^r = 
— max[— 1, max/^gar\/2 ''^b'-]- These equations can be 
solved by iteration, and their solution allows to compute 
£* (xi, t) = miuxa ^(xi, X2, t), which is found to be 

£i(xi) — > max[— 1, max hi2r]+ (2) 
i2eL2 

— N max[— 1, max Uijr] + / ma.x[0, hi^^ + ui^r] 

rert. (l2r)<£E:l2eL2 

where fi(xi) is the energy contribution of Li nodes and 
is constant relative to X2. 

A difficulty can arise if the solution is not unique and 
different solutions have different energies. In this case, 
only one of the solutions will correspond to the actual 
minimum. The following argument (confirmed by numer- 
ical investigations) suggests that this is not a problem. 

The MS equations are closed for cavity fields with sup- 
port in {—1, 1}, and also for cavity fields with support in 
{— 1, 0, 1}. Solutions with other supports can exist for fi- 
nite size instances and for appropriate initial conditions, 
but we have verified numerically that they disappear in 
the infinite size limit, so we shall ignore them. 

Let us consider (as in [T7| for the non-bipartite case) 
the uniform ensemble of instances with poissonian de- 
gree distribution and average degree c, in the infinite 
size limit. The average fraction of cavity fields ui^r 
that take the value -1-1 satisfies the equation p" = 
exp[— cexp(— cp" )]. Also (1— p!i),p+ and (l—p'l) must 
satisfy the same equation. In the case of bipartite graphs 
p^jl can be different from (and p'l from ), and this 
is a notable difference relative to the non-bipartite case. 
In any case, p^j: and p'^i are determined from p"^ and . 

For c < e the equation admits an unique solution, 
which implies that = (1 — ) and p'^ = (1 — p'i). 



meaning that the cavity fields have support over {—1, 1}. 
This unique distribution of cavity fields will correspond 
to an essentially unique fixed point of MS: it is possible 
that some disconnected components (with finite size) ad- 
mit several fixed points, but the fixed point of the 0{N) 
component (which dominates the energy) is unique. This 
is confirmed by numerical simulations. 

For c > e, the situation is more complicated: the equa- 
tion X — exp[— cexp(— cx)] admits 3 solutions, and p!f. 
can be different from 1 — (and p!f. from 1 — p^). The 
condition p" -|- p!i < 1 implies that p" < 1 — p" , so the 
total number of solutions will be at most 6. Only some of 
these possible solutions will correspond to positive values 
of the energy (and the other ones can be dismissed), and 
only one of them will be the correct one. We have stud- 
ied in detail the case for c = 5: the number of solutions 
corresponding to positive energy is 3, and remarkably the 
value of the energy is the same for all of them. One of 
the solutions has support on {—1,0,1} , corresponding 
to the 1-RSB case in the non bipartite case |T7|, and the 
remaining two have support in {—1, 1}. 

On finite size instances, we have verified numerically 
that these 3 fixed points can always be obtained by chos- 
ing appropriate initial conditions. Their energies are 
close to each other, but not exactly the same, and the 
correct one is always the largest. 

We conclude from this discussion that the energy com- 
puted from ([2]) is correct for instances extracted with 
poissonian degree distributions with c < e and approx- 
imately correct for instances with c > e. It must be 
noted, however, that the reduced instance to be solved 
in the second stage is not necessarily poissonian, as the 
probability that a node in R is matched to a node in Li 
can be correlated to its degree. Moreover, it is possible 
that some small disconnected components have multiple 
solutions, that combined with the 3 solutions of the gi- 
ant component give a larger number of fixed points, but 
these will always have energies that are approximately 
equal. We shall neglect these possible issues, comforted 
by our numerical results. 

In the following we shall give the explicit computations 
for the case where the support of u and h is {—1, 1}, but 
not for the case where the support is {—1,0,1}: even 
though we have implemented both cases, we have veri- 
fied that the energies of the solutions obtained are almost 
exactly the same for all connectivities; however, the ex- 
pressions for case {—1, 0, 1} are much more complicated, 
and the running times are much longer. 

Computing the average relative to t. To proceed with 
the computation of the average in ([T]), we note that the 
energy ([2| is a sum of local terms over Xi, h and u, 
so that its average can be computed with a procedure 
similar to Survey Propagation (SP) [TTHI3] . This cor- 
responds to a simplified Belief Propagation (BP) for the 
variables Xi, h, u and t, where Xi, h and u are subject to 
hard constraints implementing the matching conditions 
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and the MS update equations, and where t are subject 
to external fields that force them to take the marginal 

Pit) ^Ui,eL,^[th = ^] =Ui,eL,Ph- 

We introduce the probabilities Uir — Pt['Uir = 1] prop- 
agating from left to right, and Hir = Pt[/i/r — 1] propa- 
gating from right to left. The SP-like equations for Uir 
and Hir are: 

Uir ^Plllil Hir') , Hir = H Uvr) (3) 

r'^r I'^l 

Equations ([s]) can be derived by observing that 
Uir = P [t; = 1] P [-max(-l, maxr'^r ^/r) = 1| = 1], 
and similarly for Hir- 

The average minimum energy is then computed by 
averaging ^ over t using Uir = Pt[uir = 1], and 
Hir^-Ptihir = 1]: 



£:*(xi) =Etr(xi,t) = ^ 



Pi 



2Y[il-Hir)-l 



211(1 -Uir) 



2Y,HlrUlr (4) 



(/r) 



where the term £i(xi) is included and represented with 
the convention Uw = Hir = ^ir VZ e Li. For example, 
the contribution from a vertex Z g L2 is if the vertex is 
present and if all the incoming values of hir are —1, which 
happens with probability p; Y[r(^~ Hir)', the same contri- 
bution will be — 1 if the vertex is present and if there is at 
least one incoming value of hir equal to -1-1, which hap- 
pens with probability pi [1 — YI^{1 — Hir)]. The average 
of the contribution is then pi [2Y\^{1 — Hir) — !]■ The 
average of all remaining terms is computed similarly. 

Notice that the "naive" application of BP to the prob- 
lem defined over the variables Xi, h, u and t (subject 
to the appropriate external fields), in which one would 
consider the pairs (hir, uir) as single joint variables, with 
cavity probabilities Pi_j.r[(/izrj ""ir)] and Fr^ilihir^uir)], 
would lead to the wrong result for c > e. In fact, when 
the number of fixed points of the MS equations depends 
on t, the naive procedure would give to each of them a 
weight proportional to p{t) while the correct weight is 
p{t)/nt, where is the number of fixed points corre- 
sponding to a given t and p{t) is its probability. This is 
achieved with the SP procedure we introduced. 

Minimizing relative to Xi. We can then proceed to 
minimize this energy, using again MS. We consider the 
messages Uir and Hir as variables of a new problem and 
introduce the cavity messages Uir{U) = logP[[/;,r = U] 
propagating from left to right and Hir (H) = log P [Hir — 
H] from right to left. Notice that if Z G ii, we will have 
Uir = xir G {0, 1} satisfying the matching constraints 
J2r Uir < 1, while if ^ G L2 we will have Uir & [0, 1] satis- 
fying the SP update equations ([3]), and similarly for Hir- 
The continuous distributions over messages associated to 
X2 variables can be discretized for numerical purposes. 
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FIG. 1. Average energy density vs. average connectivity of L 
nodes. The three lines correspond (from top to bottom) to a 
greedy algorithm (assign xi as if t = 0) , to a "smart" greedy 
algorithm (find the maximum- weight matching on the full in- 
stance with weights on the nodes equal to their probabiUty 
to be available and assign xi accordingly), to the SP-derived 
algorithm (with h and u with support on {—1,1}), and to the 
offline lower bound of the optimum (with prior knowledge of 
t). The vertical line is at c = e. Each point is an average 
of 50 to 100 instances, with error bars smaller than the point 
sizes. The instances have \Li\ = 1000 and IL2I = \R\ = 2000, 
with pi distributed uniformly in ]0, 1]. 



equations for the 
obtained as usual 



messages Uir 
for MS, i.e. 



The update 
and Hir are 

Ulr{Ulr) = max -Eir{Ulr,Hir') + Y,r'^r'^lr'{Hir') 

where Eir{Uir, Hir') is the sum of the terms in ([4|) 
containing Uir, and where the maximisation is over 
the values of the incoming messages {Hir' r' r} 
subject to the appropriate constraints; the update of 
Hir{Hir) is obtained similarly. We don't report these 
equations for brevity. All these maximizations can be 
performed efficiently by exploiting their associativity. In 
order to improve the convergence of the algorithm we 
also introduce a reinforcement term for the messages 
associated to edges (hr) with li G Li [TBI I19j . 

These equations can be solved by iteration starting 
with uniform initial conditions. These are the only mes- 
sage passing equations that need to be solved numeri- 
cally. At the fixed point, the values of Uir and Hir pro- 
vide the optimal values of Xi by setting xir = 1 if and 
on\yii[Uir{l)-Uir{0)] + [Hir{l)-Hir{0)]+2 > 0. Once 
xi has been assigned and the realization of t has been 
extracted it is easy to perform the minimization over X2 . 

Numerical results and comparison with other meth- 
ods. Figure [1] shows some results obtained with the SP- 
derived algorithm in the case where h and u have sup- 
port on {—1,1}. The case where h and u have support 
on {—1,0, 1} gives results that are very close to these. 

In order to give quantitative evidence of the potential- 
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pending on p). For c = 3.5 and above, the time scaling 
of CPLEX worsens significantly: for p = 10, the running 
time increases dramatically with |Li |, and for |ii | = 1000 
CPLEX was not able to attain an optimum under a cut- 
off of 24 hours even for p = 2. Note that the SP-derived 
algorithm employs around one minute, (ii) The method 
was also successfully applied to a stochastic version of 
the maximum weight independent set problem |20j , when 
the node's contribution to the total weight is uncertain 
but its distribution is known. This is a relevant prob- 
lem in communication networks with some interference 
constraints [H]. Details will be given in [5^ . 
RZ acknowledges the EU grant n. 265496. 
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FIG. 2. Each point represents an instance with IL2I = \R\ = 
2|I/i| and pi distributed uniformly in ]0, 1]. Top Left: Energy 
vs. number of samples p obtained by CPLEX for c = 2.5 
and |Li| = 1 000, compared with the one computed by the 
SP-derived algorithm. The energies and error bars were com- 
puted by resampling over 10 000 samples. Bottom Left: 
CPLEX time (seconds) vs. number of samples in the same 
instances. Top Right: CPLEX time as a function of |Li 
for several values of c and p = 10. Bottom Right: Best fit of 
CPLEX times in Bottom Left with f{x) = bx°- gives a ~ 2.35. 



ities of our approach for real world problems, we have 
made two final studies: on the one hand we have com- 
pared the performance with state-of-the art method, and 
on the other we have applied the method to problems 
which are NP-hard even in the deterministic setting. In 
both cases, the results perfectly corroborate our expec- 
tations, (i) We compared the SP-derived algorithm with 
two other standard approaches. The first is a greedy 
strategy solving a weighted matching based on p{t): 
even though it is very fast, its solutions are much worse 
(Figure [T]). The second is called stochastic program- 
ming. It consists in extracting p realizations t^, . . . , f ^ 
p{t) and then solving miuxi X]f=i n^ii^x' '^(^i; 1') — 

minxj^xi,....x^ SiLi '^(^17^21 1*) using OR techniques like 
linear relaxations complemented with branch-and-bound. 
Note that this minimization problem is NP-Complete[Tl]. 
We employed two well known tools for this task: iLog 
CPLEX, a commercial, industrial strenght linear/integer 
programming software from IBM, and lp_solve, an open 
source alternative. Although qualitatively similar, re- 
sults with lp_solve were uniformly worse than the ones 
of CPLEX, so we will not report them. We observe that 
the results depend strongly on p and on the average de- 
gree c. As expected, for fixed c the quality of the solution 
improves as p increases, but the running time becomes 
larger. For c up to around 2.5, CPLEX seems to be 
able to solve the problem in polynomial time in both p 
and iV, but either it is much slower than the SP-derived 
algorithm or it gives a significantly higher energy (de- 
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